/*  How to use fft routines

Forwardfft()

Input to Forwardfft() is through the arguments xr[][] and xi[][] 
of dimension SIZE x SIZE.  The local working variable x[]
is used of dimension 2 x SIZE.  Ouput is from the 
variables xr[][] which is the real part and xi[][] which
is the imaginary part.  Each of these variables is of dimension 
SIZE x SIZE where SIZE is the size of the largest square image you wish to work 
with (e.g. float gxr[SIZE][SIZE]).
gmaxx and gmaxy are global variables indicating the size of the input targa image where gmaxx < SIZE and gmaxy < SIZE.  
To take the forward dft, just put gpicval into xr[][] and zero's into xi[][], then call Forwardfft().

Backwardfft()

Input to Backwardfft() is through the arguments xr[][] which is the real part,  and xi[][] which is the imaginary part.  
Each are dimensioned as SIZE x SIZE (e.g. float xr[SIZE][SIZE]). The local working variable x[] is used of dimension 2 x SIZE.  
Ouput is from the arguments xr[][] which is the real part and xi[][] which is the imginary part. 
Each of these variables is of dimension SIZE x SIZE where SIZE is the size of the largest square image you wish to work with.

gmaxx and gmaxy are global variables indicating the size of the input targa image where gmaxx < SIZE and gmaxy < SIZE 

*/


void Forwardfft(float xr[SIZE][SIZE],float xi[SIZE][SIZE] )
{
        float scale,x[2*SIZE];
        int i,j,k,l;
        scale=gmaxx*gmaxy;

        printf("forwardfft: \n");

        for (k=0;k<gmaxy;k++) 
		{
                for(j=0;j<gmaxx;j++)
			{
                        x[2*j]=xr[j][k];
                        x[2*j+1]=0.0;
                	}
                Fft(x,gmaxx,1.0);
                for(j=0;j<gmaxx;j++) 
			{
                        xr[j][k]=x[2*j]; 
                        xi[j][k]=x[2*j+1];
                	}
        	}
        printf("horizontal pass fft done\n");

        for(j=0;j<gmaxx;j++)
		{
                for (k=0;k<gmaxy;k++) 
			{
                        x[2*k]  =xr[j][k];
                        x[2*k+1]=xi[j][k];
                	}
                Fft(x,gmaxy,1.0);
                for(k=0;k<gmaxy;k++) 
			{
                        xr[j][k]=x[2*k]/scale;
                        xi[j][k]=x[2*k+1]/scale;
                	}
        	}
        printf("vertical pass fft done\n");
        printf("forwardfft() done\n");
}



/*___________________________________________________________*/           

void Backwardfft(float xr[SIZE][SIZE],float xi[SIZE][SIZE])
{
        float scale,x[2*SIZE];
        int i,j,k,l;
        scale=gmaxx*gmaxy;

        printf("backwardfft: \n");

        for (k=0;k<gmaxy;k++) 
		{
                for(j=0;j<gmaxx;j++)
			{
                        x[2*j]=xr[j][k];
                        x[2*j+1]=xi[j][k];
                	}
                Fft(x,gmaxx, -1.0);
                for(j=0;j<gmaxx;j++) 
			{
                        xr[j][k]=x[2*j];      
                        xi[j][k]=x[2*j+1];
                	}
        	}
        printf("horizontal pass inverse fft done\n");

        for(j=0;j<gmaxx;j++)
		{
                for (k=0;k<gmaxy;k++) 
			{
                        x[2*k]  =xr[j][k];
                        x[2*k+1]=xi[j][k];
                	}
               	Fft(x,gmaxy, -1.0);
                for(k=0;k<gmaxy;k++) 
			{
                        xr[j][k]=x[2*k]/scale;
                        xi[j][k]=x[2*k+1]/scale;
                	}
        	}
        printf("vertical pass inverse fft done\n");
        printf("forwardfft() done\n");
}


/*___________________________________________________________*/

void Fft(float data[], int num, float dir)
/*      Discrete fft algorithm
 *      forward transform called by Fft(data,num,1.0);
 *      inverse transform called by Fft(data,num,-1.0);
*/
{
        int array_size, bits, ind, j, j1;
        int i, i1, u, step, inc;
        float sine[513],cose[513];
        float wcos, wsin, tr, ti, temp;
        bits = log(num)/log(2) + .5;
        for (i=0; i<num+1;i++)  
		{
                sine[i]=dir*sin(3.141592654*i/num);
                cose[i]=cos(3.141592654*i/num);
        	}
        array_size=num<<1;
        for(i=0;i<num;i++)  
		{
                ind=0;
                for (j=0;j<bits;j++)  
			{
                        u=1<<j; ind=(ind<<1)+((u&i)>>j);
                	}
                ind=ind<<1;j=i<<1;
                if(j<ind)  
			{
                        temp=data[j]; 
			data[j]=data[ind];
                        data[ind]=temp;
			temp=data[j+1];
                        data[j+1]=data[ind+1];
                        data[ind+1]=temp;
                	}
        	}
        for(inc=2;inc<array_size;inc=step)  
		{
                step=inc<<1;
                for(u=0;u<inc;u+=2)  
			{
                        ind=((long)u<<bits)/inc;
                        wcos=cose[ind];wsin=sine[ind];
                        for(i=u;i<array_size;i+=step)  
				{
                                j=i+inc;j1=j+1;i1=i+1;
                                tr=wcos*data[j]-wsin*data[j1];
                                ti=wcos*data[j1]+wsin*data[j];
                                data[j]=data[i]-tr;
                                data[i]=data[i]+tr;
                                data[j1]=data[i1]-ti;
                                data[i1]=data[i1]+ti;
                        	}
                	}
        	}
}
/* _________________________________________________________________ */